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1 . SUMMARY 

The computer code for calculating web temperature distribution 
has been expanded to provide a graphics output for (aT)" in addition to 
numerical and punch card output. The new code was used to examine 
various modifications of the J419 configuration and, on the basis of the 
results, a new growth geometry was designed. Additionally, several 
mathematically defined temperature profiles were evaluated for the 
effects of the free boundary (growth front) on the thermal stress 
generation. 

Experimental growth runs were made with modified J419 configura- 
tions to complement the modeling work. A modified J435 configurauion 
was evaluated. 


complete 


description of the computer code is included as an appendix. 
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2 . INTRODUCTION 


Silicon dendritic web is a single-crystal silicon ribbon 
material which provides substantial advantages Tor low-cost manufacture 
of solar cells. A significant feature of the process is the growth from 
a melt of silicon without constraining dies, resulting in on oriented 
single-crystal ribbon having excellent surface features. In common v/ith 
other more classical processes such as Czochralski growth, impurity 
rejection into the melt permits the use of less pure "solar grade" 
starting material without significantly affecting cell performance. A 
unique property of the dendritic web process is the growth of long 
ribbons of controllable width and thickness which not only facilitates 
automation of subsequent processing into solar cells, but also results 
in high material utilization since cutting and polishing are not 
required. 

On the present contract, three broad areas of work are 
emphasized: 

1. The development of thermal stress models in order to 
understand the detailed parameters which generate buckling 
stresses. The model can then be used to guide the design of 
improved low-stress web growth configurations for 
experimental testing. 

2. lixperiments to increase our understanding of tlie efCects of 
various parameters on the web growth process. 

3. The construction of an experimental web growth machine which 
contains in a single unit all the mechanical and electronic 
features developed previously so that experiments can be 
carried out under tightly controlled conditions. 

Thus, the principal objective of this work has been to expand 
our knowledge and understanding of both the theoretical and experimental 
aspects of the web growth process to provide a solid base for 
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substantial Imptovements in both area throughput and web crystal 
quality. 

During this reporting period, the new thermal model was used to 
examine various modifications of the J419 configuration. These results 
combined with experimental evaluation formed the basis for a new growth 
design. Work continued to optimize the J435 configuration for steady- 
state growth and a modified slot geometry was evaluated. 
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3. TLCliNiCAL PROGRESS 


3.1 Modeling 

3.1.1 Model Development 

During the previous quarter, a new model was developed to 
calculate the temperature profile in a growing web crystal. The 
objective of the new model was to provide a amch higher resolution in 
the representation of the lid and shield geometry of the growth system 
which had been treated as more or less lumped blocks in the previous 
version. As discussed in the previous report, the new model provides 
for an exact geometrical representation of each lid, shield, and spacer 
of a growth configuration. In this report, we discuss in more detail 
the input parameters available for the model as well as the output 
options. 

Input options that are available for the model fall into two 
catagories: selection of the operating options of the model such as 

type of output, number of cases, mode of integration, etc. and, second, 
the parameters of the configuration being analyzed. A description of 
the individual data cards as well as a complete listing of the program 
are included as an appendix to this report. 

Output options available for the program fail into three 
categories: 1) numerical output, 2) graphics output, and 3) punched 

card output-: The first and third output types were previously 

available; the second output was added when it became apparent that the 
second derivative of the thermal expansion coefficient times the 
temperature was a useful parafiieter for assessing the effect of changes 
in the shield elements on the thermal stress. 
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An example of the numerical output is shown in figure 1. The 
first column lists the position on the web (the growth front is at 0 
cm). The second column gives the calculated tempfiratuce. The third and 
fourth columns are tin- first and second derivatives of the temperature, 
and the fifth column is (aT)"^ All distances are in centimeters and 
temperatures are in °K. In addition to this data, two other sets of 
data are listed at the top of the printout. The first, labeled VV, is 
the partial growth velocity of the web resulting from the heat lost from 
the web crystal itself. This differs from the total growth velocity 
which is observed in experiments by a contribution whicn results from 
some of the latent heat being dissipated to the supercooled raelt.^^^ An 
estimate of this total growth velocity is included as one of the input 
parameters to the model. 

The second set of data included is an estimate of the critical 
yield stress for the web corresponding to the first five temperature 
points. The yield stress is simply an estimate based on an empirical 
equation from the work of Graham et al.^^^ 

°YP “ 2.57 X 10“^^ exp(49459/T) Mdyn/cm^ (1) 

These critical-yield stress values can be used later when actual stress 
distributions are generated by the finite element WHCAN calculations. 
Although the viscoelastic phenomena responsible for the observed 
residual stress in the crystals is extremely complex, these data provide 
a first estimate of possible effects. The first four columns in the 
balance of the numerical data are more or less self-explanatory in being 
the position, the temperature, and its first two derivatives. The last 
column, (aT)", is of importance since it is basically the generating 
function for the thermal stress. Although the relationship between 
(uT)’’ and stress is not straight forward near the interface where the 
free boundary exerts a strong influence, it is relatively simple 
elsewhere where the Bo ley and Weiner approximatlon^'^^ applies: 
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Figure 1. Example oE numerical output from new web temperature computer 
code. 
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0^ “f (3y^ “ w2) (aT)" (2) 

where E is Young’s modulus and w is the half width of the ribbon. 

Through the relationship of T" with the heat loss from the ribbon, it is 
possible to correlate the behavior of this parameter with changes in the 
geometry of the lids and shields. 

Realization of the usefulness of the (uT)" parameter led to the 
addition of a graphics output capability to the code. Since the end 
application was to compare the (aT)'* function with the geometry being 
analyzed, the graphics output presents both a representation of the lid 
and shields and the concomit?i|i‘ t (ctT)" curve. If the geometry does not 
change (lid, shield, and interface position), then up to throe different 
(aT)" curves can be represenU<^d on the same plot. A representative 
output is shown in Figure 2j the plot of the lid and shields has been 
accented in the figure for the sake of emphasis. 

The final output option from the program is a set of punched 
cards giving the nodal temperatures for use with the WECAN stress 
calculations. In this option, there are several sub-options depending 
on whether a two-dimensional stress calculation or a three-dimensional 
buckling calculation is to be done. Further, in the two-dimensional 
calculation, there is a choice as to whether quadratic or cubic elements 
are employed. Also, some choice is available as to some of the 
geometric features of the finite clement grid. 

Thus, the program to calculate the web temperature distribution 
has a great deal of flexibility in both the lid and shield configuration 
which can be analyzed, and in the options for presenting the results of 
its calculations. The flexibility of both the input and output greatly 
enhance the ability of the program to assist in the analysis and design 
of dendritic web growth configurations. 
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Figure 2. Graphical output from new web temperature computer code. 
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3»1«2 Modeling; AppllcatlotiB 

3. 1.2.1 "Roalistic" Growth ConflBurationB 

The first application of the new web teraperatvire model was an 
evaluation of variations of the J419 configuration (iHgure 3). Of 
course> the individual lids and shields wore modeled rather than the 
lumped parameter "model” shown in the figure. For some of the 
geometries evaluated, only the (uT)" curves were generated and, in these 
cases, the initial peak of the curve and the final, small, and maximum 
were the features of greatest interest. They are related respectively 
to the extremum In the y-stress at the growth front and to the maximum 

in AOj^ (Uiffereneo in x-stress between center and edge) which occurs 

further up the crystal. In turn, .hose stresses are related to the 
reKidual stress in the web and to the tendency of the web to buckle. 

One of the first parameters studied was the effect of work coil 
height us it affected the lid and shield temperatures. Temperature data 
were obtained from thermocouple measurements of the bottom lid and the 
bottom and top shields; other shield temperatures were obtained by 
interpolation. Although the total range of temperature change in the 
lid was only 18“k, the affect of changing coil position was primarily 

near the growth front; Initial y-stress and (partial web velocity). 

The emu* omit mil shield temperature clianges were over but the AOj^ 

and the associated (uT)" peak were changed very little. These results 
indicate that the design of the lid itself is very Important in 
controlling both residual stress and growth speed, but doe« not have a 
great deal of control on the buckling. 

As a follow-up, a beveled lid slot and a straight ltd slot were 
compared. The beveled configuration was about 10% faster, but the 
initial y-stress increased about 30%. Again, minimal changes occured in 
the distant Ao,^ peak. 

One gaometrlcal parameter whicli did influence tlie Ao^ peak was 
the heiglit of the snleld stack. Increasing the separation between 
shields to give a higher stack not only moved the peak f’'i‘ther away from 
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th.‘ interface but decreased its magnitude. All other factors being 
equal, this should result in a much decreased tendency to buckle. 

Other geometrical parameters investigated were the width of the 
lid slot, the view angle through the shield stack, and the effect of 
solid spacers between the shields, etc. It was found that decreasing 
the slot width should Increase the growth speed as does opening the viev/ 
angle. All of these factors also tend to Increase the initial y-stress, 
so that there is a trade-off in arriving at any final design. 

On the basis of all the factors studied, a new proposed growth 
configuration was designed and constructed. This design, called the 
"J46U" configuration (after its first trial run), was anticipated to 
have low buckling tendency, although the residual stress was 
uncertain. In fact, as discussed elsewhere, this proved to be a very 
successful configuration. 

3. 1.2. 2 ^Synthetic Temperature Profile 

In addition to modeling existing or hypothetical lid and shield 
configurations, stresses were calculated for two "synthetic" temperature 
profiles that were defined ns purely mathematical functions. The goal 
of these test cases was to obtain an estimate, and perhaps at least a 
semiquantltatlve representation of the free-boundary end effect. The 
rep.iou near the growth front is of prime Importance in the generation of 
plastic deformation in the web crystals and yet is a region of least 
theoretical guidance in relating the temperature profile to the thermal 
stresses. 

The first synthetic temperature profile was a rerun of the 
constant stress profile discussed in a previous report. Although the 
previous results were close to those predicted from simple analytical 
considerations, there were several aspects where they differed. At 
first, this was attributed to numerical effects In the WliCAN program, 
but then it was realized that the finite element grid had edge elements 
which modeled the bounding dendrites and would cause the stress fields 
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to depart from the Ideal cade of a flat ribbon. 'X’he temperature profile 
w«*8 re-analyzed using a mesh for a simple ribbon, and the results were 
now in perfect agreement with the theoretical predictions. 

As before, the effect of the free boundary extended into the 
crystal for a distance approximately equal to the full ribbon width. As 
an approximation, the x-stress varied as 1-exp (-x^/w^), where w * 
ribbon half width. Using this relation as an admittedly approximate 
guideline, it is possible to assess the relative effects of features in 
the (uT)" profile near the interface. 

The second synthetic temperature profile assumed an exponential 
decay of the (aT)" profile from its value near the interface. This 
particular function was chosen since it is a reasonable approximation to 
the behavior observed in models of any realistic lid/shield 
configurations, although in those cases, the exponential behavior is 
found only in the first 5 mm or so, The stress distribution from the 
synthetic temperature case was in reasonable quantitative agreement with 
the stress distributions in the interface region of the more realistic 
cases, although of course the more distant features were quite 
different. In both the synthetic and realistic cases, the actual x- 
stress magnitude near the boundary was far smaller than might be 
anticipated, in agreement with the boundary effect of the constant 
stress case. These results indicate that the initial (aT)" peak is 
responsible for much of the thermal stress behavior at the interface, 
but also tliat there are other factors such as the free boundary effect 
which must be considered. It is also not yet clear to what effect the 
peak height and the characteristic decay length of the exponential are 
involved in the stress generation. These factors need some additional 
runs to clarify the behavior. The final results should be applicable to 
the design of lid configurations for faster growth with lower residual 
stress . 
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3 . 2 lixperlme»tal Web Growth 

3.2.1 Introduction 

The process of developing a functional new growth configuration 
follows a three-stage progression. The first stage is the thermal 
modeling, the results of which generate a design for a growth 
configuration. Stress and buckling models are applied as appropriate 
when the (uT)" results warrant the effort. The design is then 
fabricated into hardware and tested experimentally. The objective of 
this second phase of the progression is to experimentally verify the 
stress behavior predicted by the model, i.e., to find how wide the 
crystal can grow before deformation occurs. At this stage, long slots 
which provide mcl t temperature profiles compatible with crystal widths 
of about b.5 cm are used. Lid and shield temperature measurements are 
mode to verify consistency with the model, and modifications are made as 
required to produce the desired temperature distribution in the vertical 
direction. This stage involves a cross interaction between modeling and 
experiment. The third stage of development is the adaptation of the 
low-stress configuration for semi-automated growth. This phase is 
largely empirical, guided by experience. It involves the incorporation 
of melt replenishment and width control provisions into the design, and 
the determination of optimum sliielding and coil position for good growth 
at constant width and melt level, i.e., steady-state growth. Only the 
designs which continue to show promise through the first and second 
phases reach the third phase. 

3.2.2 Growth Experiments 

A series of variations of the JA19 configuration were tested in 
order to evaluate the effects on growth of changes in top shield spacing 
and minor variations in slot geometry and coil position, etc. Lid and 
shield temperature measurements were made to complement the modeling 
work. Growth parameters evaluated were growth velocity and residual 
stress, in addition to the width at which buckling was initiated. 
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In general, Che only modification of the basic design which had 
a positive effect on growth was a small bevel on the lid slot. Other 
variations either had little observable effect on growth behavior or 
generated a negative result. For example, opening the slots in the top 
shields increased the growth velocity but caused the crystal to 
degenerate at narrower widths. The model indicated that increasing the 
height of the top shield stack should further reduce buckling stresses, 
but this did not seem to have a significant effect with the J419 lid 
configuration. Finally, on Che basis of modeling results and the 
experimental correlations, it was concluded that a new lid slot geometry 
combined with an extended shield stack should produce significantly 
lower buckling stresses than the J419 configuration. The required 
hardware has been fabriesfced and the design will be fully evaluated 
during the next reporting period. However, preliminary Indications are 
that the new design is indeed a major step forward in stress reduction 
and width enhancement. 

In parallel with the above effort, experimental growth runs with 
the J^35 configuration (which combines the low-stress aspects of the 
JAiy with the width-limiting capabilities of the J98M3) were continued 
in the N-furnace. The objective was to optimize the furnace parameters 
for steady-state growth. The parameters to be determined Include coil 
height, end shield positions, and melt level. 

having adjusted these parameters, a number of crystals was grown 
3 to 4 meters long with the width consistently saturating at 3.3 cm, the 
design target for this configuration. However, at this width, a 
considerable portion of either edge of the web is growing in the dog 
bone region of the slot, i.e., the full face of the web does not see a 
uniform slot width. This results in a visible change in the web surface 
along tne edges as compared with tlie central region. In an effort to 
improve surface uniformity, a modified lower lid was fabricated without 
dog bones La the growth region, i.e., a straight sLot. 
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The surface of web grown with this straight slot was extremely 
flat and uniform. However, the web crystals widened slowly and were 
subject to pullout if not closely monitored by the furnace operator. In 
one run> a 4.3 meter long, 3.4 cm wide crystal was grown, but this was 
the exception rather than the rule. Thus, although the straight slot 
produces a beautiful crystal surface, it is difficult to grow a crystal 
at full width. We plan a further modification such that at full width, 
only the dendrites see the open lid environment, while the full face of 
the web sees a uniform slot. This arrangement should improve growth 
stability while maintaining uniform surface quality. 




4. CONCLUSIONS 

The new thermal model has proved a valuable tool in its ability 
to accommodate complex lid and shield configurations. Its application 
to modifications of the J419 configuration along with experimental input 
has led to a new design for a low-stress growth configuration in a 
systematic manner. 

Semi-automated web growth has been demonstrated with the J435 
configuration in the N-furnacej although we are not yet satisfied that 
this basic configuration is fully optimized for long-term steady-state 
growth. Further trimming of this configuration will continue in 
parallel with evaluation of the new design aimed at much greater crystal 
width. 


ORIGINAL 

OF POOR QUALITi' 


5. PLANS AND FUTURE WORK 


The computer code will be used for further evaluation of a new 
lid design and the bucRling behavior of the configuration will be 
calculated. The new lid design will be evaluated experimentally and 
modified as required. 


6 . NEW TECHNOLOGY 


No new technology is reportable for the period covered. 
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9 . PROGRAM COSTS 


Costs 


Costs 


27,682 Previous $1,232,077 

3,415 This Quarter 162,088 

31,097 Cumulative $1,394,165 
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appendix I 

Web TeraperaCure Computer Model 
1. INTRODUCTION 


Th-srmal stress and buckling modeling of dendritic web crystals 
requires temperature distribution data along the growing web* Because 
of the difficulty of measuring web temperature, we Instead compute it 
from the furnace lid and shield geometry and temperature distribution. 
This computation is performed by a computer code called "RIBBON" which 
integrates the required heat transfer equation for the given radiative 
web environment. A description of this program and data input follows. 

2. CONSTRUCTION OF MODEL 


2 * I Differential Equation 

The purpose of the RIBBON program is to determine the 
temperature, T, of the growing silicon web as a function of the 
distance, x, from the lower edge of the furnace lid. RIBBON 
accomplishes this task by integrating the heat conduction equation: 


p CpV 


dT , d la dT) 
dx " dx 4 dx-* 



where p = density 

Cp » specific heat 
V = web pull velocity 
a - 318 W/cin 
b = web thickness 

q « heat flux from one side of the web. 


lA-1] 


O A 


One of the major tasks of the RIBBON program Is to calculate the 
geometric form factors for the terra "q". A detailed discussion of this 
computation is found in the last quarterly report <D0E/JPi;.“955843/82/6) . 
For the purpose of recognizing the degree of nonlinearity of equation 
A-1, we note that q can be expressed in the form 

q « oET*^ d- f(x), 

where o is the Stephan-Boltzman, E is the eraissivity of silicon web, and 
f(x) is a function of position x on the web and the geometry and 
temperatures of the lid and shields of the furnace. Above 20 cm from 
the lid, the web enters a chimney and thus no longer "sees" the lid and 
shields. For a simple approximation, we say Chat it "sees" only the 
ambient temperature, fg* Thus, for x r 20 cm,, 

q , oe (T^ - [A-2] 

The longer the ribbon grows, the closer its temperature approaches the 
ambient. This fact together with the initial condition that the web 
starts to grow at the silicon melt temperature, T^ “ 1685°K, gives us 
the boundary conditions 


T(x„) » T„ 

'o' m 

T(“>) « T^ . [A-3] 

The position x = x^ is the growth front and has a value equal to the 
negative of the input parameter LIN in the program. 

Differential equations with boundary conditions are more 
difficult to solve than those with initial conditions; they generally 
must be solved Iteratively. Since most numerical Integration routines 
are written for first order equations, we transform the second order 
equation (A-1) into two first order equations with the substitutions: 
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U2 - p Cp V T^/a 
32 « eo T^^/(ab) 

TT(1) « T/T^ 

TT(2) - II (1/T) [A-4] 

Equation A-1 can be expressed now as a system of first order equations: 

TTP(l) “ TT(1) • TT(2) 

TTP(2) » U2 TTP(l) + $2^ 

where 

TTP(i) - 4^ [TT(i)] i * 1,2 

Cl A 

and 

Q » 2q/eo . [A-5] 

The nonlinearity of these equations makes them unstable as small 
errors in the steps of the integration are quickly magnified in the 
succeeding steps. We tried several different methods of numerical 
integration but found the simple fourth order Runge-Kutta method to be 
the most stable. Even so it is necessary to use -double precision to 
obtain reasonable results. In practice, an initial slope of the 
temperature is guessed. If the slope leads to a curve which gives below 
the ambient temperature, then the slope is increased for the next 
guess. If the resulting curve goes above the silicon melt temperature, 
then the initial slope is decreased. This iteration continues until the 
double precision accuracy of the choice of initial slope is exhausted; 
in other words, there is no way to choose an initial slope between the 
high slope and the low slope since they are identical to 16 places. 
Generally, the integration curve does not blow up (or down) until 
20 - 30 cm. In this case the integration from 0 to 10 cm (the length of 
the buckling finite element model) is fairly accurate. If it blows up 
before this length, this may no longer be true. (A smaller integration 


22 


ORlGSNAi PAQS 
OF POOR QUAirrY 


8tep size, UO, may improve this problem.) Also, if the integration 
remains between the melt and ambient temperature for a longer length, 
the integration may not be accurate because other values of the initial 
slope may lead to different temperature curves which also reraoin 
bounded. The boundary condition, T («®) ■ T^, must then be applied. 

Thus, it becomes necessary to examine the asymptotic expansion of 
equations A-1 and A-2. While the numerical integration of equation A-1 
or equation A-4 starting with initial values at the melt interface will 
is be accurate for small values of x, the asymptotic expansion can be 

expected to be accurate for large values of x. Hopefully, their regions 
of accuracy will overlap; if this is rot the case for some problems, 
then some approximation to an intermediate solution might be required. 

2.2 Asymptotic Expansion 

For large x, it is most convenient to change the differential 
equations A-1 and A-2 into a system of first order equations with the 
following substitutions} 


Thus , 


a •= p C[ " Tq/u 
Y * 2eo Tg^^/aba^ 
yi = T/Ta 


[A-61 




dx ““^1^2 

"^'^2 4 

iT ‘ y I y-i 


lA-71 


If one attempts a formal power series solution of equation A-7 in 
negative powers of x (so that they are finite as x approaches infinity), 
one would obtain the trivial solution = 1 and Y 2 “ This solution, 
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howcvoi*, dQO» noU filvo Uhe Bomci’al nayiopbotic CKpitnalon of equittilon 
A-7 . Theoi'y (Woll'Brtng Waoon, AaympUotLc Expmiolona Cor Ordinary 
Dtf foroniiai Hquatlnno, intorHcionci* > N.Y., 1965) ohowo Chat the gcnoral 
aolutlon ,l8 a Cunctlon of uwo parumoUera: G exp ax(l“K M>16 y)/ 2 and 

1) C/’p ax(l‘l’/i‘iloT)/'-i> where C and 1) are arbitv’ary couutnnin* Blnco the 
aeeond parameUor bocomea Inflnluo mi x upproachoo InClivity, we let 
I) vanifih and look for funckiono of the first parameter. The aubstlkution 

5 u exp rnxa-/?PT6Y)/a] 


trnnafurms the equation A“6 Into 


y^^ i; (1 ~ iWlTy) » 2 y^y^ 

y^'5 U - /r-PISY) = 2 y,^ + 2Y (y/' - D [A-ll] 


where the prliiioa represent dirferentiatlon with reopcct to the variable 
f, . We can now seek a power aeriefi aolntion of equation A-8 in the, form 


^1 " ‘^10 ‘^11 '’12 

2 

y.^ “ a.2| (Gfi) 5* (^C) d* ... 

where “» 1 from the boundary condition and n^ may alao be choaen 
equal to unity ainao G to an arbitrary eonauanu. Equating the 
GoeCficlentB of (GO*^ on both aideB of equation A-8, we obtain 


n-1 


n a (1 ~ /fPTSY) “2 ). a a 


In 


mwO 


Im 2(n-m) 


tA-101 


n-1 „ (n/2) -1 

n a, (1 - /T+foT) ^ I a, , + 2y [), , ^ +2 ] b li ) 

2n li» 2(a-nO n/2 m n-m' 
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fo? n even or 




n-1 (n-n/2 

^ Jo “1" Jo 


[A-10] 


for n odd, 


where 

\ = ^(1/2)' + J, \i "Ui-J) - ^ "'''^■' 

(l-l)/2 

’ Jo i 


b = 1 
0 

From these equations, the aj^^^ and &2\x determined from the values of 
aj^^ and a 2 „i, where m < n. In this way, we find the coefficients of the 
asymptotic e'. vans ion A-9. The radii of convergence of these power series 
may be found from the formulas: 

R = Lini 1 an/^n+l ' exists 

and Lim sup 1 a^^ j^/n = a « 1/R 

n -»■ 00 

The next steps in the procedure are to choose a value of C and then a 
value of ? > 20 cm. such that equation A-9 converges. From equations A-4, 
A-6, and A-9, initial conditions for equation A-5 may be found at the 
value of X corresponding to that of Equation A— 5 can now be integrated 
backward to x = x^ using the same Runge-Kutta method as before. These 
steps are interated by increasing or decreasing the choice of C 
accordingly as T(Xq) is less tlian or greater than T^. 


25 


3. COMPUTER PROGRAM INPUT/OUTPUT OPTIONS 


The input parameters of the RIBBON program are divided into two 
sets — the first to define the geometry and the second everything else. 

In this way, Calcorap plots may be obtained for any number of geometric 
configurations in each run. In each plot, up to three different graphs 
(black, green, and red) of the second derivative of aT can be obtained for 
different nongeoraetric parameters of web growth run.* A sample Calcorap 
plot is illustrated in Figure 2. Besides the oT second derivative curves, 
it diagrams the lid and shield geometry. 

The following input records are read in the first data set. The 
data are read in free field format with either spaces or commas separating 
them. Variables beginning with I-N are integer, (except for LIN) and all 
others are double precision. 

RECORD 1; 

IGRAPH = Input 0 if there are to be no Calcorap graphs in the run, 
otherwise input 1. 


RECORD 2; 

NS = Number of geometric elements. These include the lid, the 
shields, and their separating gaps. 

JQC = Parameter to correspond to the type of WECAN element 
= 2 for quadratic elements 
= 3 for cubic elements 

JN = Number of WECAN elements in the x direction (generally has been 
set to 23). 

NE = Number of WECAN elements in the y direction (generally has been 5 
for two-dimensional elements and 7 for three-dimensional ones). 


*Here a stands for the thermal expansion of silicon and is represented 
by ALPO + ALPI * TEMP in the program; it is unrelated to the a in 
equations A6 and A7 . 
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NG = One plus the number of data sets for a given geometry. 

UN = Distance (cm) of the growth front below the lower edge of the 

lid. 

EMAG = Magnification of the JNth element in the x direction relative to 
the first WECAN element at the growth front (generally has been 
8 ). 

ELL = Length of ribbon modeled on WECAN in centimeters (generally has 
been 10). 

RECORD 3: 

11(1), I * 1, -NS = Height (cm) of ith horizontal surface above the lower 
surface of the lid. 

RECORD 4: 

Y(I), 1=0, ..., NS+1 = Half width (cm) of the central gap of the ith 

horizontal surface. 

The H and Y parameters are illustrated in Figure 3. This completes the 
first data set. 

The second data set consists of the following records: 

RECORD 1: 

IS =1 for forward integration 

= -1 for backward integration from asymptotic expansion (not yet 
implemented) . 

HO = integration step size (generally set to 0.01). 

KOUT = 6 for printed temperature data in WECAN stress input formalism 
= 7 for punched WECAN data. 

KOU = 0 for no plot of the uT second derivative. 

= L for a plot 

EPS = eraissivity of silicon web. 
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A 

B 

V 

FOMIN 
POMAX ^ 
CMIN 

CMAX 

RECORD ; 
TS(I),I 


chemal conductivity of silicon multiplied by its temperature 
(318 W/cm). 

web thickness (cm). 

web pull velocity (cm/min) 

minimum initial value of TT(2). If not known, set to -1. 

maximum initial value of TT(2). If not known, set to -0.01. 

minimum value of the parameter in equation (A-9). If not known, 
try 1. 

maximum value of the parameter C. If not known, try 20. 

= 1, ... , NS+1 = Temperature (degrees Kelvin) of the ith 
geometric element. 


These RECORDS 1 and 2 may be repeated up to a total of three sets. The 
data sets 1 and 2 can be repeated an indefinite number of times. 
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APPENDIX II 
Computer Code Listing 


2RUN, /RNPT JSS05,09FA0REA0Y20,SCHRUBEN, 5, 100/1500 
aHDG RIBBON 


aPTNis 
FTN 8R 


5?sr 


08/25/82-20:50 

1. C SILICON WEB - ANALYSIS OF HEAT LOSS OF NOVING WEB 

I: 585fe^i‘;«|^yitoVE!S;^{Sg,^§78!?55lHco=^o^,»co:^o),»s<o=1oo>. 

%IMENSI0N*TT<2) ,TS4 (20) ,T( 100) ,TP(100) ,TPP(100> »ATPP(100) ,Y2(( 
INTEGER “^IXSCAO) iISR<0:60i20) 

INTEGER TITLE(11),LAB£LXC2)»LABELY(2),IC0L0R(4) 

REAL XR(150),YR(A,150) jXLEN 
THE ABOVE ARE GRAPH PARAMETERS 

COMMON /CNI/ U2,BETA2,ISRtTS4,Y2,H»ICHIH,NS,NS1 


4. 

1: 

7. 

8 . 

9. 

10 . 

11 . 

12 . 

13. 

14. 

15. 

18 = 

19. 

20 . 
21 . 


24. 

25. 

ik 

28. 

29. 

30. 

31. 

32. 

33. 

34. 

35. 

36. 

37. 

38. 

39. 

40. 
41 . 
42. 
43i» 

44. 

45. 
4 6. 

47. 

48. 

49. 

i?; 

52. 

53. 

54. 

55. 

56. 

57. 

58. 

59. 

60. 
61. 
62. 

63. 

64. 

65. 

66 . 

67. 

68 . 



ATIV't'E 

00*^/ 


OF' 


ALP','HA*T‘ 


DATA 1C0L0R/I,1|2t3/ 

Kl= 5 
K0= 6 

IGRAPH = 0 FOR NO CALCOMP 

= 1 FOR SOME CALCOMP „ 

NG IS ONE PLUS NUMBER OF DATA SETS FOR A^GIVEN GEOMETRY (4 MA) 
NS IS NUMBER OF SHIELDS PLUS LID SECTIONS PLUS GAPS INBETWEEN 
TS(I) IS THE TEMPERATURE OF THE ITH SHIELD OR GAP ^ 

H(I) IS THE HEIGHT OF THE ITH HORIZONTAL SURFACE ABOVE THE LOl 
SURFACE OF THE LID 

YU) IS THE DISTANCE OF SHIELD EDGE TO WEB 

(YCO) IS HALF-WIDTH OF MELT SURFACE, Y(1> IS HALF GAP OF LlD,i 
YCNS) IS DISTANCE OF INNER EDGE OF„TOP^SHIELD TO WEB. 

Y(NS+1) IS DISTANCE OF OUTER EDGE OF TOP SHIELD TO WEB) 

JQC=2 FOR QUADRATIC ELEMENTS 

=3 FOR CUBIC ELEMENTS • - 

IS'-I FOR FORWARD INTEGRATION 
=-1 FOR BACKWARD INTEGRATION 
JN IS NUMBER OF ELEMENTS IN A ROW ALONG THE WEB AXIS 
LINgU O^SJANCE OF GROWTH FRONT BELOW LID 

"^^EMAG^is’^MAGNIFICATION OF LAST ELEMENT RELATIVE TO FIRST NEXT 1 
ELL IS LENGTH OF WEB TO BE MODELED 
WE ASSUME ELL < 20 

■ FOR 2D ELEMENTS 
FOR 3D ELEMENTS 
1685. DO 


ASSUME 
NE= 5 


TS 

TM= 1685.00 
TMR= 1./TM 
LAMB= 1804. DO 
TM4= TM**4 
CP= 0.981100 
RHO= 2.30D0 
SIGMA= 5.67D-12 
H<0)= 0 
X(0)= -1 
1SR(0,Q)= 0 
X0= X(0) 

ISR(0,1)= 1 
XS(0)= XO. 
REAOCKI,*) IGRAPH 
IF (IGRAPH .GT. 0) 
YR(1,1)= -1 
XR(1)= -1 
XR(2)= -1 

DO 998 ISUPER=1,10 
READ<K1,*,EN0=999) 
WRITE(KQ,5) 

5 FORMAT! 

1 

NS1= NS+1 
H(NS1)= 20 


CALL PL0TIN(1.,1 .) 


NS,JQC,JN,NE,NG,LIN,EMAG,ELL 
" LIN = ‘',E12.5,' EMAG = *',E12.5,' pLL=‘',E12.5) 


> 14 , 
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RIBBON 

1 

1 

1 

1 

1 

1 

1 

2 

2 

1 

1 

1 

1 

1 

1 

2 

2 

2 

2 

2 

2 

3 

4 
4 
4 


69. 

70. 

71. 

72. 

73. 

74. 

75. 

76. 

77. 

78. 

79. 

80. 
81. 
82. 

83. 

84. 

85. 

86 . 
87- 
88 . 

89. 

90. 

91. 

92. 

93. 

94. 

95. 

96. 

97. 

98. 

99. 
100 . 
101 . 
102 . 

103. 

104. 

105. 

108. 

109. 

IIS: 

112 . 

113. 

114. 

115. 

116. 

117. 

118. 

119. 

120 . 
121 . 
122 . 

123. 

124. 

125. 

126. 

127. 

128. 

129. 

130. 

131. 

132. 

133. 

134. 

135. 

136. 

137. 

138. 

139. 

140. 

141. 

142. 

143. 

144. 


READIKI,*)^ 
WRITE (KQ, 7) 

formatPoh 


^ CH{I),I = 1,NS). 


IjNS, 

.#.vci*«.5)) 

READtKI,*) CYa>,I = 0,NSl>, 

WRITE<I<Q.8) CyCl) ,I=0,NS1 ) 

8 FORMAT(*OY= ^,/,(D14.5)> 

00 10 1= 0.NS1 
Y2(I)= YU>**2 
10 CONTINUE 

YR(1,2)» SNGL(YCO)) 

ALPO= 2.8457E~6 
AUP1= 9.796E~1Q 
ALP2= 2+ALP1 
K« 0 

DO 50 1= 2, NS 
11 = 1-1 
XLS= XO 

XLS IS THE LOWER BOUND OF TS(I) RADIATION ON WEB 
HI= H(l) 

Yl= Y(l) 

DO 20 J= 1,11 

IF <YI .GT. V<J)> THEN 

TEMP= (H(J)AYI-Y (J)*HI) FCYI-Y(J)) 

IF (TEMP .GT. XLS) XLS= TEMP 

END IF 
20 CONTINUE 

IF >(XLS .LE. XO) THEN 
1SRCQ,1)= 1 

ELSE 

ISR(0,I)= -1 
DO 25 J=K,0,-1 
J1= J+1 

IF (XLS .GE. XS(J)) THEN 
XS(J1)= XLS 
IXS(J1)= I 
GO TO 26 

ELSE 

XS(J1)= XSCJ) 

IXS(J1)= IXS(J) 

END IF 

25 CONTINUE 

26 K= K+1 
END IF 

XUS= 20 

XUS IS THE UPPER BOUND OF TS(I1) RADIATION ON WEB 
HI= H(I-2) 

YI= Y(I1) 

DO 30 J= I, NS 
J1= J-1 

IF (YI .GT, V(J)) THEN 

TEHP= (Y1*H(41)-HI*Y(J))/(YI-Y(J) ) 

IF (TEMP .LT. XUS) XUS= TEMP 

END IF 
30 CONTINUE 

IF (XUS .LT. 20) THEN 
DO 35 J= K,0,-1 
J1= J+1 

IF (XUS .GE. XS(J)) THEN 
XS(J1)= XUS 
IXS(J1)= II 
GO TO 36 


ELSE 


XS(J1)= XS(J) 
IXS(J1)= IXS(J) 


35 

36 

50 


60 


END IF 
CONTINUE 
K= K + 1 

END IF 
CONTINUE 
IB= 0 
1 = 0 

DO 70 L= 1,K 

IF (XS(L)-X(I) .LT. 1.D-4) THEN 
IXSL= IXS(L) 

ISR(I,IXSL)= -ISR(I,IXSL) 
60 TO 70 

30 




RIBBON 




145. 

146. 

147. 

148. 

149. 

150. 

U5: 

153. 

154. 

it?: 

157. 

158. 

159. 

160. 
161. 

ftl: 

164. 

165. 

166. 

167. 

168. 

169. 

170. 

171. 

172. 
17J. 

174. 

175. 

176. 
177., 

178. 

179. 
ISO. 
181. 
182. 
183. 

Itt: 

186. 

187. 

188. 

189. 

190. 

191. 

192. 

193. 

lit: 

196. 

197. 

198. 

199. 

200 . 
201 . 


2( 

204. 

205. 

206. 

207. 

208. 

209. 

210 . 
211 . 
212 . 

213. 

214. 

215. 

216. 

217. 

218. 

219. 

220 . 


end if 

I1» I 
1 = 1+1 

_ DO 65 J» O.NS, ^ ^ 

65 1SR(I.J)= ISR(I1,J) , 

IF (H(13)-XS(L) .GT. 1.D-4) THEN 
^ X(1)=.XSCL) 

IXSL= IXSCL) 

ISR(1,IXSL)= -1SR(I,IXSL) 

EJ-SE . . . . 

XCl)s HUB) 

IB= IB+1 
1SR(I,0)» IB 
GO TO 60 

END IF 
70 CONTINUE 

DO 80 L= IB, NS , 

IFCHCL) .LT. X<I)> GO TO 81 

1 = 1 + 1 
X(I)= K(L) 

DO 75 J= 1,NS 
75 ISR(I,J)= ISR(I1,J) 

80 ISRC1,0)= L+1 

81 CONTINUE 
LAS= I 

IF (XCI) .NE. 20) THEN 


LAS= 1+1 
X(l ■ ‘ 


:las)= 20 

END IF 
LAS= LAS+1 
X(LAS)= 30 

WRITE CKO, 1000) (X( J) ,J-0 .LAS) 
WRITE(K0,1001)(XS(J),IXSCJ),J=1,K) 

WRITE C KO, 1 002 ) CC1SR(L,J),J=0 ,NS) ,L = 0,I) 

DO 90 1= 1,LAS 
K= I 

1F((XCI)+UN).GT. 1.D-8) 60 TO 91 

90 CONTINUE 

91 CONTINUE 
K0= K 

RA= EMA6 + *(1.D0 /(JN-D) 

EL= EUL*(1~RA)/(1-RA**JN) 

JNN= jac*JN+i 
KK= 2 

DO 110 KKI= 0,K 
KKKS= KKI+1 
KKK- KKI 

IF(CH(KKI)+LIN) .GE. 0. ) 60 TO 111 

KK= KK+1 

XRCKK)= SNGLCHCKKI)). 

YRC1,KK)= SNGLCYCKKI)) 

KK= KK+1 

XR(KK)= SNGLCHCKKD) 

YR(1,KK)= SN6LCYCKKKS) ) 

110 CONTINUE 

111 KKS= KK 
KKS1= KK+1 

DO 997 IC= 2,NG „ ^ 

READ (K I,*) IS,HO,KOUT,KOU,EPS,A,B,V,POMIN,POMAX,CMIN,CM 
WRITE (KO, 82) I S ,HO,KOUT ,KOU , EPS , A ,B , V ,P0MIN ,POMAX ,CHIN , 
KOUT = 6 FOR PRINTED STRESS DATA 
= 7 FOR PUNCHED STRESS DATA 
KOU > 0 FOR GRAPH 

82 FORHATC'O IS = '’,I3,' H0=',E12.5, ^ , 

1 ' K0UT = ',13,' K0U=‘',I3,' EPS= ',E1 2.5,'^^A= ',E1 

2 E12.5,' V= ',E12.5,/, " P0MIN= ',D25.18,' P0MAX= 

3 , 't ,CMINf ',025.18,* CHAX“ ',D25.18> 

REAOCKI,*) (TSCI),I=1,NS1) 

WRITE(K0,6) (TSCI) .1=1 ,NS1> 

6 FORMAT('OTS=',/, (014.4)) 

DO 83 1= 0.NS1 

TS(I)= TS(I)*TMR 
TS4(I)= TS(I)**4 

83 CONTINUE 

IS1= M0D(1-IS,JQC) 
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RIBBON 


221 . 


224. 

Ill: 

227. 

228. 
229. 

111 ! 

Ill: 

234. 

235. 

236. 

237. 

238. 

239. 

240. 

241. 

242. 

243. 

244. 

245. 

246. 

247. 

248. 

249. 

250. 

251. 

252. 

253. 

lit: 

256. 

257. 

lit: 

260. 

261. 

262. 

263. 

264. 

265. 

266. 

267. 

268. 

269. 

270. 

271. 

272. 

273. 

274. 

275. 

lit: 

278. 

279. 

280. 
281. 
282. 

283. 

284. 

285. 

Ill: 

288. 

289. 

290. 

291. 

292. 

293. 

I?t: 

296. 


400 

409 


410 


400 


.NE. 1) 


GO TO 102 


-1) X2= XC1) 


411 

412 


HO* H0*IS 
TA* TS(NSI) 

V* V/6D 

BETA2* EPS*SIGMA*TH4/(A*a) 

BETAS 0SQRTC2*BETA2J 
U2* RHO*CP*V*TM/A 
DO 101 J» 1,500 
IFLAGs 0 
K* KO 

OX* EL/JQC 

IF (IS .EQ. -1) GO TO 
TT(D* 1 
XX* -LIN 

POOO* .500*(P0MAX+P0MIN) 

IF (POOO .EQ. POO .AND. J 

???I.E“?8o 

KK» KKS+1 
KKK= KKKS-1 
11 * 1 
12= JNN-1 
X2= -LIN 
ICHIM* 1 
GO TO 409 
CONTINUE 
continue 
TEMP* TT(1)*TM 
TEMPP* TT(2)*T£MP 
DO 450 1 = 11 ,12 ,IS 
X1= X2 
X2= X1+0X 

IF (I .EQ. 2 .AND. IS .EQ. 

KFLA6* 0 
T(I)= TEMP 
TP(I)= TEMPP 
XS(I)= XI 

IF (IS*CX(K)-X2)) 410,411,412 
CALL FCNJ(X1,X(j() ,TT,TTP,HQ,K) 

IF (IFLAG .NE. 0) 60 TO 600 
XR(KK)= SNGL(XI) 

TEMPPP* TEMP*TTP(2)+TENPP*TEMPP/TEMP 
ATEM* (ALPp+ALP2*TEMP)*TEHPPP+ALP2*TEMPP*TEMPP 
TEMP* TT(1)*TM 
TEMPP* TT(2)*TEMP 

YRCIC.KK)* AMAXU-1 ..AMIN1 (8 . ,SNGL (ATEM ) *500 . ) ) 

IF (KFLAG .EQ. 0) THEN 

TPP(I)* TEMPPP 

ATPP(I)= ATEM 

END IF, 

IF (H(KKK) .EQ. XI) THEN 
KKK* KKK+1 
KK= KK+1 

XR(KK)= SNGL(XI) 

YR(IC,KK)= VR^1C,KK-1) 

END IF 
KK* KK+1 
KFLAG* 1 
XI* X(K) 

K* K*IS 

IF (IS*(X(K)-X2)) 410,411,412 
CALL FCNJ (XI ,X2,TT,TTP ,HO,K) 

IF (IFLAG .NE. 0) GO TO 600 

TEMPPP* TEMP*TTPC2)-*'TEMPP*TEMPP/TEMP 

ATEM* (ALP0+ALP2*TEMP)*TEMPPP+ALP2*TEMPP*TEMPP 

TEMP* TT(1)*TM 

TEMPP* TT(2)*TEMP 

XR(KK)* SNGLCX1) 

YR(IC,KK)= AMAXK-1 .,AMIN1 (8 . ,SN6L (ATEM )*500 . ) > 
IF (H(KKK> .EQ. XI ) THEN 
KKK* KKK+1 
KK* KK+1 

XR(KK)= SH6L(X1) 

YR(IC,KK)= YR(IC,KK-1) 

END IF ^ 

KK* KK+1 

IF (KFLAG .NE. 0) GO TO 450 
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III: 

5?S: 

301. 


30S. 


i; 


309. 

II?; 

312. 

313. 

si 

316. 

317. 
316. 

319. 

320. 

321. 

322. 

323. 

T 3 i. , 

ft. ■ 

325. 

326. 

327. 

328. 

329. 

330. 

331. 

332. 

333. 

334. 

335. 

Ill: 

338. 

339. 

It?: 

342. 

343. 

344. 

345. 

346. 

347. 

348. 

349. 

II?: 

352. 

353. 

354. 

355. 

356. 

357. 

Ill; 

360. 

361. 

362. 

363. 

364. 

365. 

366. 
367a 
368. 

IJ?; 

2 / 1 . 

372. 


TPPCl)« TEMPPP 
ATPP<I)» ATEM 

450 IF (MOOCI.JQC) .EQ. 2S1> DX» 0X*RA**IS 
CALL FCNl(X2,TT»TTP,K)f 
13= I2+1S 
TCI3)« TT(1>#TM 

lsCl)= X2 
TP(I)= TT<2)*T<1> 

ATPpJ j“=^!ALp5+ALP2*T?in*TPP<I)+ALP2*TP(I)**2 

xr(kkJ* SNGLCX2) 

YR(IC.KK)= AMAXK-1 .,AM1N1 C8 . »SNGL <ATPP Cl ) ) %500 . ) ) 

IF CIS .EQ. -1) GO TO 600 
DO 500 JK= K.LAS 
IF (JK .EQ. LAS) ICHIM= 0 
CALL FCNJ(X2,XCJK),TT.TTP,H0,JK) 

I|^(I^LA^ .NE. 0) GO TO 600 

CONTINUE 
CONTINUE 

WRlTE(KO,10g2) IFLAG 
WRITE<KO,1005) XI ,X2,TTf P00,C 
FORMat(3d25.18) 

IF CIS .EQ. -1) GO TO 95 
IF (IFLAG) 93.92,94 
GO TO 102 
P0HIN= POO 
GO TO 101 
94 POMAX= POO 
GO TO 101 
CONTINUE 
CONTINUE 

VV= -60*A* POO/CRHO*LAMB) 

WRITE(KO,7100> VV ^ 

7100 F0RI*IATC'1',SX,'VV= ',E14,7> 

DO 849 r=i;S 

Y1ELD= 2.57E-11*EXPC49459/TC1)) 

URITE(KQ,8490) YIELD 
F0R{JAT(" Y1ELD= ^,E15.8) 

DO 852 1=1, JNN 

WR1TECK0,1000)XS(I),TC1) ,TPCI),TPPCI),ATPP(1) 

DO 850 1=1 ,JN 
11= JQC*(I-1)+1 
T1= TCII) 

T2= TClI+1) 

51: 55J!tl! 

DO 800 J=1,NE 

11= J+CI-1) ANE 

IF fJQC .EQ, 3) GO TO 795 

WRITE CKOUTy 1006) II.TI ,T1 ,T3,T3,T1 

IF (NE aEQ. 5) THEN 

WR1TECK0UT,1008) T2,T3,T2 

£ In S E 

URirE(KOUT,1007) T1 , T3 ,T3 , T1 ,T2 ,T3 ,T2 ,T1 ,T1 ,T3 
WRITECK0UT,1008) T3,T1 ,T2,T3,T2 

END IF 
GO TO 000 

795 WR1TE(K0UT,1006) 1 1 ,T1 ,T1 , T4 ,T4 , T1 
WRITECKOUT, 1007) T2 ,T4 ,T3,T1 ,T3 

aoo 

850 CONTINUE 

YRCIC.1)= 8 
DO 360 KKI= 2,KKS 
YR(IC,XK1)= YRC1C,XKS+1) 

860 CONTINUE 
997 CONTINUE 
ISW= -1 
K= KKKS-1 

DO 865 KKI= KKSI.KK 
IF (XR(KKI) .NE. SNGLCHCK))) THEN 
YR(1,KKI)= YR(1,KK1~1) 

ELSE 

K= K+ HAXuCO.ISW) 
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RIBBON 


173. 

!?§: 


yR<1tKKI)« SN6L<V(K>) 
ISW» -ISW 

END IF 

865 CONTINUE 

DO 866 1= 1,KK 

866 WRITECKO,1009) . XRCn.tV 

XLEN» l.-SNGLCLIN) +ELL 
IF (KOU .6T 


CYRCJ,I),Jc1,NG> 


1 


998 CONTINUE 

999 CONTINUE 
IF (I6RAPH .6T 
F0RMAT<8D15«8) 
FORMAT(7(D15.8iI 3)> 
FORf-ATOOU) 
FORNATn6,T2X,5F12.4, 

luuf F0RMAT(18X,5F12.4,' ♦ 

1008 FORMATC18X,5F12,A) 

1009 FORMAT(8E15.8) 

CALL EXIT 

END 


0) CALL GRAPH<XR,YR,NG,4,KK,XLEN.9.»0,0,TITLEf 
LAB£LX«6(LASELY»6|0fOrICOLOR) 


1000 

]001 

1002 

1006 

1007 


0) CALL PLOT(0.,0.»999) 


+ ') 


V 
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RIDDON 


END IF 
11* I, 

1 * r + 1 

DO 65 a* OiNS^ ^ ^ 

65 ISR(I,J)» ISR<I1,J> 

IF (H(I3)-XS(L) .GT» 1,0-6) THEN 

usi. miC)”''-’ 

ISRU,IXSL>« -ISR<I,IXSL) 

ELSE 

ISR(I,0>» IB 
GO TO 60 

END IF 
70 CONTINUE 

00 80 L* IB, NS 

IF <H(L) .LT. X<D) 60 TO 81 
II* I 
I» 1+1 
X(I)« HCL) 

00 75 Jf l.NS, 

75 ISR(I,J)» ISRdliJ) 

80 ISR(1,0)» L+1 

81 CONTINUE 
LAS* I 

IF (XilJ ,Ni. 2 ( 3 } THEN 


X^LAsK^ZO 

END IF 
LAS» LAS+1 

WRI^E^Ko5?000) (X(J) ,J»0, ,LAS) 

WRJ,TJ:(K0,10gi) (XS(J),IXS(J>,J*1,K) 
WR|iC(K0,1002)(ClSR(L,a),a«0,NS),L»0,I) 

DO yd I» 1.LAS 
X* I 

IF(<X<I)+LIN).6T. 1.0-8) 60 TO 91 

90 CONTINUE 

91 CONTINUE 
K0« K 

RA= EMAG**(1.D0 /(JN-1)> 

EL* ELL*<1-RA)/C1-RA**JN)( 

JNN= JQC+JN+1 
KK« 2 

DO 110 KKI= 0,K 
KKKS* KKl+1 
KKK® KKI 

1F((H(KKI) + HN) .6E. 0. ) 60 TO 111 

KK® KK+1 

XR(KK)* SNGLCH(KKD). 

YRd.KK)* SNGL(YIKKI)) 

KK= KK+1 

XR(KX)= SNGLCHCKKD) 

YRII.KK)* SNGLCYCKKKS)) 

110 CONTINUE 

111 KXS= KX 
XKSis KK+1 

DO 997 IC= 2,NG 

READ(KI, + ) lS,MOi>KOUT.KOU,EPS.A,B,V,POniN,POMAX,CMIN.CNAX 
WRITE (KO,82) I S .HO.KOUT ,KOU , EPS , A ,B ,V , PONIN ,PONAX ,CriIN , CH AX 
KOUT = 6 FOR PRINTED STRESS DATA 
= 7 FOR PUNCHED STRESS DATA 
KOU > O^FOR GRAPH ^ 

82^F0RHAT( Q | *5 , , ^ g-, 2,5 ^ ',E12.5,‘ 

REA0(KI,*J CTSCI),I®1,NS1) 

WRITE(K0,6) CTSCI),1=1,NS1) 

6 FORMAK'OTSs'./.IdU.A)) 

DO 83 1= 0.NS1 

TS(I)= TS(I)+TMR 
TS6CI)= TSCI)**A 
83 CONTINUE 

IS1= M0D(1-IS,JQC) 


rJ^O «LAS) 
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